In-class Exercise 6

In this exercise we learn how to work with Spatial Autocorrelation

1. Setup

pacman::p_load(sf, sfdep, tmap, tidyverse, DT)

hunan <- st_read(dsn = "data/geospatial", layer = "Hunan")

Reading layer `Hunan' from data source 
  `/Users/jeffery/Projects/Y4S1/IS415/S0methingSimple/IS415-GAA/in-class_ex/ex06/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 88 features and 7 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 108.7831 ymin: 24.6342 xmax: 114.2544 ymax: 30.12812
Geodetic CRS:  WGS 84

hunan2012 <- read_csv("data/aspatial/Hunan_2012.csv")

hunan <- left_join(hunan,hunan2012) %>%
  select(1:3,7,15,16,31,32)

# nb: A neighbor list object created by st_neighbor
# style: "W" for standardized weights
# .before = 1: Insert into the front
# Returns a table with the neighbors that can be viewed

wm_q <- hunan %>%
  mutate(nb = st_contiguity(geometry),
         wt = st_weights(nb, style = "W"), .before = 1)

datatable(wm_q)

moranI <- global_moran(wm_q$GDPPC, wm_q$nb, wm_q$wt)

glimpse(moranI)

List of 2
 $ I: num 0.301
 $ K: num 7.64

K is average neighbours found

# Performs a basic test
global_moran_test(wm_q$GDPPC, wm_q$nb, wm_q$wt)

    Moran I test under randomisation

data:  x  
weights: listw    

Moran I statistic standard deviate = 4.7351, p-value = 1.095e-06
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.300749970      -0.011494253       0.004348351 

p-value: < 0.05, enough statistical evidence we can reject null hypothesis, does not conform to random distribution (95% sure) - Fail to reject null hypothesis if greater than 0.05 no point proceeding Moran I: Greater than 0 suggest clustering in the data, but is a relatively low clustering

# Always best practice to set seed before simulation
set.seed(1234)

Global Moran I Permutation

global_moran_perm(wm_q$GDPPC, wm_q$nb, wm_q$wt, nsim = 99)

    Monte-Carlo simulation of Moran I

data:  x 
weights: listw  
number of simulations + 1: 100 

statistic = 0.30075, observed rank = 100, p-value < 2.2e-16
alternative hypothesis: two.sided

p-value: Is smaller Moran I: Is the same

Local Moran I Permutation

lisa <- wm_q %>%
  mutate(local_moran = local_moran(GDPPC, nb, wt, nsim = 99), .before = 1) %>%
  unnest(local_moran) # Turn it into a table form

datatable(lisa)

ii: local moran i p_ii: p-value with base method p_ii_sim: Based on simulation mean: label hotspots median: Use if there is significant no. of positive or negative skew Plot out skewness to see if you want to use median

tmap_mode("plot")
pii_m <- tm_shape(lisa) +
  tm_fill("p_ii") + 
  tm_borders(alpha=0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Local Moran's p-value", main.title.size = 1)

ii_m <- tm_shape(lisa) +
  tm_fill("ii") + 
  tm_borders(alpha=0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Local Moran's I", main.title.size = 1)

tmap_arrange(pii_m, ii_m, ncol = 2)

lisa_sig <- lisa %>%
  filter(p_ii < 0.05) 

tm_shape(lisa) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
  tm_shape(lisa_sig) +
  tm_fill("median") + 
  tm_borders(alpha = 0.4)

Hot Spot Cold Spot Analysis

# MUST Use distance inverse distance, Further are smaller
wm_idw <- hunan %>%
  mutate(nb = st_contiguity(geometry),
         wts = st_inverse_distance(nb, geometry, scale = 1, alpha = 1),
         .before = 1)

HCSA <- wm_idw %>%
  mutate(local_Gi = local_gstar_perm(GDPPC, nb, wts, nsim = 99), .before = 1) %>%
  unnest(local_Gi) # Turn it into a table form

datatable(HCSA)

gi_star value:

LISA use cluster and outlier G*: Hot and cold spot

h_pii_m <- tm_shape(HCSA) +
  tm_fill("p_sim") + 
  tm_borders(alpha=0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Local G* p-value", main.title.size = 1)

h_ii_m <- tm_shape(HCSA) +
  tm_fill("gi_star") + 
  tm_borders(alpha=0.5) +
  tm_view(set.zoom.limits = c(6,8)) +
  tm_layout(main.title = "Local G* I", main.title.size = 1)

tmap_arrange(h_pii_m, h_ii_m, ncol = 2)

HCSA_sig <- HCSA %>%
  filter(p_sim < 0.05) 

tm_shape(HCSA) +
  tm_polygons() +
  tm_borders(alpha = 0.5) +
  tm_shape(HCSA_sig) +
  tm_fill("gi_star") + 
  tm_borders(alpha = 0.4)